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Abstract. Many experiments, and in particular gravitational wave detectors, 
produce continuous streams of data whose frequency representations contain discrete, 
relatively narrowband coherent features at high amplitude. We discuss the application 
of digital Fourier transforms (DFTs) to characterization of these features, hereafter 
frequently referred to as lines. Application of DFTs to continuously produced time 
domain data are achieved through an algorithm [7], hereafter referred to as EFC§, 
for efficient time-domain determination of the Fourier coefficients of a data set. We 
first define EFC and discuss parameters relating to the algorithm that determine its 
properties and action on the data. In gravitational wave interferometers, these lines 
are commonly due to parasitic sources of coherent background interference coupling 
into the instrument. Using GEO 600 data, we next demonstrate that time domain 
subtraction of lines can proceed without detrimental effects cither on features at 
frequencies separated from that of the subtracted line, or on features at the frequency 
of the line but having different stationarity properties. 

PACS numbers: 04.80.Nn, 02.30.Nw, 07.05.Kf 
1. Introduction 

The purpose of this paper is to describe and validate a method for characterizing and 
removing coherent backgrounds, hereafter referred to as lines, in continuously acquired 
data. There has been much work in this area, including several studies of techniques 
intended specifically for application to data from gravitational wave detectors. See for 
example references [1, 2, 3, 4, 5]. These existing methods make very limited use of 
Fourier coefficients, because the usual methods for obtaining these coefficients through 
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the fast Fourier transform (FFT) algorithm [6] run into practical difficulties where it is 
required that the Fourier coefficients be equally accurate estimators of line parameters 
for all data samples, and that the estimate of line parameters be sufficiently fast to keep 
pace with the data sampling rate. 

In this paper, we describe the application of an alternative to the fast Fourier 
transform, referred to throughout the paper as EFC, for obtaining the Fourier coefficients 
continuously updated at the data sampling rate. This algorithm was described in 
1966 by Halberstein [7], and is the subject of U.S. and E.U. patents. See [8, 9], for 
example. The authors exploit this algorithm for the monitoring of coherent backgrounds 
in gravitational wave detector data. We demonstrate that EFC is significantly faster 
than real time running on contemporary computers, and equally accurate for all data 
samples. We continue by providing a proof of principle of EFC using data from a ground- 
based gravitational wave detector, and by studying the response of the algorithm to 
injected signals. 

A key application of a robust real time capable line removal tool is the generation of 
triggers for networks of interferometric gravitational wave detectors in real time for the 
purposes of providing data for follow up of triggers with optical or radio observations. 
Applications of EFC in this area are more widespread than the removal of lines. A 
continuously maintained monitor of magnitude and phase of coherent signals in the 
data has many applications in gravitational wave interferometery, including calibration 
of the detectors, monitoring interferometer subsystems such as high Q suspensions, and 
studies of bicoherence between lines and features at other frequencies in the same data. 

2. Line parameterization with Fourier transforms 

By a line we mean an oscillating feature in the data whose amplitude Ai at the time, 
t i: of the i th sample can be written 



where Aq, /, and are the underlying amplitude, frequency, and phase of the 
spectral feature, which are functions of time t. The algorithm we will describe requires 
that these parameters do not vary significantly on timescales less than some cutoff r. We 
will also confine the discussion to lines having f(t) combined to a well defined domain 
for all t, and to data sampled regularly, meaning that the sampling time, t s = t n+1 — t n , 
is independent of n. 

If the frequency / is known precisely, it only remains to evaluate the amplitude A 
and phase at each time t n . A natural starting point would be the k th coefficient of 
the digital Fourier transform [14] of N successive data samples Xj centered on the n th 
data sample. 



A n — 



A cos(27r ft n + 0), 



(1) 



N-1 




(2) 



3=0 
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where k ~ Nt s f and Nt s < r. In the case where / is not precisely known, or 
corresponds to non-integer k, the values of several Fourier coefficients would be needed. 
Fourier coefficients can be extracted from the DFT of the data, which is frequently 
calculated using a fast Fourier transform (FFT) algorithm [14, 15]. The number of 
floating point operations to calculate the DFT of iV successive data points is of order 
N\og 2 N. 

This method is beset with difficulties in the case where one requires a continuous, 
sample-by-sample estimator of A, f and 0. Because fast Fourier transform algorithms 
require a block of time series data, it is not obvious how to proceed from an estimate 
of A n and <p n at the n th data sample to an estimate of the A n+ i and <p ri+ \. It is too 
computationally intensive to perform an FFT for each time sample - iV samples would 
require order of N 2 log 2 N operations. The alternative usually considered is performing 
a new FFT every m samples, where m < N. The problem in this case is the occurrence 
of discontinuities in the estimator at each new transform. 

3. The evolving Fourier coefficient (EFC) algorithm 

In the case under discussion, only a restricted number of Fourier coefficients are 
of interest, suggesting an alternative approach. One starts by generating Fourier 
coefficients of iV data points using Equation 2. It is shown in Appendix A that the 
k th Fourier coefficient, J-k{x\, X2, • • • ,xn), of time samples x±, X2, • • • ,xn is related to 
FkixQ, x 1 ,---,x N _ 1 ) by 

Fk(xi, xn) = (Tk(x Q , • • • , £jv-i) + (x N - x )) . (3) 

Using this result, the computation of the Fourier coefficient of the overlapping time 
series is achieved without the need for an additional Fourier transform. The process 
can be repeated to obtain the time evolution of the Fourier coefficient for all subsequent 
time samples. Figure 1 illustrates the algorithm for the case where the sum is over four 
successive data samples. 

We recall some properties of discrete Fourier transforms whose choice will be 
important in practical application of this technique to line monitoring and removal. 
First, the Fourier coefficient determines the line parameters over a data set of duration 
r = Nt s . This timescale should be shorter than the typical timescale for variations in 
the line parameters A, f and 0, but longer than the typical timescale for features in the 
data that you would not wish to influence the line parameterization. The duration r 
also determines the resolution bandwidth B = 1/r, the bandwidth in frequency space 
sampled by a single Fourier coefficient. If the lines under study have a width or range of 
wander Af, then the number of Fourier coefficients to be calculated is rAf = Nt s Af 
per new time sample. 

EFC requires an initial value for each Fourier coefficient of the data. We note, 
however, that if we use zeros for the values of x subtracted in Equation 3 for the first 
N timesamples, the output of the algorithm exactly equals the Fourier component after 
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Figure 1. A block diagram illustration of the operation of the EFC algorithm for 
the simple case where each FFT is over four data samples. In the example iteration 
illustrated here, the sum of the new data sample and minus a data sample taken four 
sampling times previously are added to the previous Fourier coefficient, and the result 
of this sum is multiplied by exp(i<5), where 8 = 2irk/N is the phase difference between 
successive time samples at the frequency of the Fourier component. The newest data 
sample then replaces the fifth most recent data sample in the ring buffer. 



N iterations. Hence there is no need for initial values of the Fourier components, since 
the iteration algorithm generates them. Equation 3 implies that the algorithm must 
evaluate the difference between the most recent data value and the data value x 
acquired N + 1 sampling times previously. After startup of the algorithm, the latter 
sample is unavailable until N sampling times have elapsed, hence there is an unavoidable 
start-up period of Nt s to initialize the filter during which the algorithm output is not 
useful for line characterization. The response of the algorithm to transients is discussed 
in Section 4.6. A natural and fast way to implement access to the N + l most recent data 
samples is to store N previously acquired samples in a ring buffer, where at each iteration 
the oldest stored sample Xq is read, applied in the algorithm, and then overwritten with 
the most recent sample xn- A counter is then incremented, modulo N, so that the next 
element in the ring buffer, containing the input data value to be used as x in the next 
iteration of Equation 3, is referenced. 



4. Tests on GEO 600 data 

4-1. Characteristics of the test data 

The immediate application for which the algorithm was conceived is the removal of 
coherent noise backgrounds from data produced by gravitational wave interferometers, 
such as the GEO [16], LIGO [11], VIRGO [12] and TAMA [13] instruments. For a 
general review of the science of direct gravitational wave detection the reader is referred 
to [10]. 
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Figure 2. Amplitude spectral density (ASD) of noise from the GEO 600 michelson 
photodiode channel Gl : LSC_MIC_VIS. The upper plot shows amplitude per ^/Hz as a 
function of frequency, where each bin has a bandwidth of 7.8 mHz. The dashed and 
solid curves are of data before and after the whitening prefilter discussed in Section 
4.1. The lower plot is the cumulative ASD of the whitened data, obtained by adding 
in quadrature successive bins from the bold curve in the upper plot and normalizing 
to the maximum value. Lines contributing significantly to the time domain amplitude 
of the signal are shown as large, discrete steps in the cumulative ASD. 

The tests of EFC described in this paper were carried out on data acquired 
from the GEO 600 gravitational wave interferometer through the LSC- Virgo scientific 
collabration. The GEO600 detector has been described in detail elsewhere [16]. For 
the purposes of testing this algorithm, an arbitrary 1000 seconds of data from the 
Gl :LSC_MIC_VIS channel, sampled at 16384 Hz, and acquired when the instrument was 
in maintenance mode, was used. This signal represents a measurement of the DC light 
power reflected from the Michelson at the input port, commonly referred to as the 
common-mode visibility. 

Figure 2 shows the amplitude spectral density (ASD) of the data before filtering. 
Plot (a) shows power spectral density over 7.8 mHz bins from 1 Hz to 8192 Hz, the 
Nyquist frequency of the data acquisition system. The dotted curve on this plot shows 
the ASD of the raw data. Broadband noise below 50 Hz and above 1 kHz dominates the 
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Table 1. Performance of the subtraction algorithm for the fifteen lines treated in 
the test. At each frequency (column 1), a number of bins (column 3) corresponding 
to a bandwidth (column 2) about the line centre were treated with the subtraction 
algorithm. After subtraction, the percentage shown in column 4 of the line amplitude 
before subtraction remained. Line amplitudes were determined using the frequency 
bin having the maximum ASD before subtraction. 



Frequency (Hz) Bandwidth (Hz) Number of bins Remnant (%) 
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cumulative amplitude spectral density. These components are suppressed using an IIR 
whitening filter [14] consisting of two zeros at Hz, two real poles at 50 Hz and a single 
real pole at 1 kHz. The solid curve is an ASD of the data after this preprocessing step. 

Plot (b) from figure 2 shows a cumulative ASD, obtained by squaring the ASD of 
the prefiltered data, integrating over frequency, taking the square root of the result, 
and dividing by the contents of the highest frequency bin. Line features in the data 
making the largest contributions to the overall signal amplitude appear as steps in the 
cumulative ASD. Table 1 column 1 shows the frequencies of some of these lines. 

4-2. Line removal 

The practical application of equation 1 to remove lines in the data is to estimate 
the amplitude and phase of oscillatory components of the preprocessed data at the 
frequencies of identified lines, to synthesize an oscillation having this amplitude and 
phase, and to subtract this oscillation from the data. The first step on receiving a new 
sample of data Xk is to use equation 1 to generate an updated Fourier component y^. 
The real part of this Fourier component is proportional to the wave of frequency k/Nt s 
contributing to the data. Note, however that the phase of y& represents the phase at 
the beginning of the timeseries, whereas for subtraction we require the phase half way 
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through the timeseries ||. Therefore y& is modified by a phase factor representing the 
phase shift over N/2 — 1 samples, and then the real part of the result is taken. The 
correction signal is denoted and is given by 



where N and have the same meaning as in Section 3 and 3? denotes the real 
part. One final important detail is the time lag between the data and the correction. 

is the Fourier component of input data Xk-N+i • • • Xk, therefore it best represents the 
amplitude and phase of the corresponding frequency component at data sample k — N/2. 
Therefore the actual subtraction of the synthesized wave is performed on a timesample 
N/2 sample times lagging behind the k th input sample. In this sense the line subtraction 
filter is acausal; it requires data both before and after the point of interest to determine 
the correction at that point. For post-processing removal of lines, for example as part 
of an off-line analysis chain, this is not a problem. The line subtracted data set {r^} is 
given by 



for all k. Note that the above equation subtracts a synthesized wave representing 
only one Fourier component from the data set. The method can be applied to multiple 
components, and to lines spanning more than one bin in Fourier space. The bandwidth 
of adjacent frequency bins is the reciprocal of the time duration of the equivalent Fourier 
transform. 

4-3. Computational overhead of the processing algorithm 

A practical line removal algorithm must be able to process the data at a rate at least 
as high as the production rate for the data set. Fast data channels from existing 
gravitational wave detectors sample at t s = 1/16384 s. The algorithm was implemented 
in the programming language C, and benchmarked on a MacBook Pro having a clock 
speed of 2.33 GHz and an Intel Core 2 Duo chipset, using the ANSI C time library 
and the computer's on-board clock. The time per incoming data sample per Fourier 
component analyzed was 1/1700 of t S) so the code can process 1700 Fourier components 
in parallel on a single CPU and match the incoming data rate. Because each Fourier 
component is treated separately, the algorithm could be run on a cluster in which case 
the number of frequencies treated is limited only by the number of cluster nodes. Tests 
on an intel based PC running linux yielded 1730 Fourier components in t s . 

4-4- Performance on the test data 

Performance on the test data is summarized in Figures 3, 4, and 5. Figure 3 shows the 
cumulative ASD before and after subtraction on the same scale. The r.m.s. amplitude 

| In practice the time series has an even number of samples, so we take the subtraction to occur at 
sample k-N/2, which is the leftmost of the two samples in the centre of the time series 
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Figure 3. Cumulative amplitude spectral density (ASD) of noise from the GEO 600 
michclson photodiodc channel Gl :LSC_MIC_VIS before and after subtraction of the lines 
listed in Table 1 . Both curves are scaled to the maximum of the cumulative amplitude 
spectral density of the unsubtracted data. At the Nyquist frequency, the amplitude of 
the line subtracted data is 0.54 that of the unsubtracted data. 



of the prewhitened data is reduced to 0.54 of its unprocessed value after line removal. 
Figure 4 shows the ASD of the input data with the ASD of the correction superimposed; 
as expected the correction signal is dominated by sharp peaks at the frequencies of each 
of the corrected lines. Plot (b) of this figure shows the ASD of the data before and 
after line subtraction overlaid. Plot (c) shows the ratio of the bin-by-bin ASDs of the 
data after to before line removal; each line is removed to a maximum of 24% of its pre- 
treatment value, with the largest lines (having the highest SNR in the power spectrum) 
attenuated to 0.6% of their amplitude in the input data. The 4 th column of Table 1 
shows the percentage of the pre-subtraction ASD remaining at the Fourier component 
dominating the line after subtraction. These numbers correspond to the depth of the 
dips at the line frequencies in Figure 3, plot (c). 

4-5. The effect of a noise background 

Noise in the subtracted signal enters through the noise background in the frequency 
bin being monitored. This noise background is considered here as broadband Johnson 
noise having a flat amplitude spectral density across the bin. As discussed in Section 
3, the bandwidth of a single bin is given by B = 1/r where r is the duration of the 
timeseries used to estimate each Fourier coefficient. If the amplitude of the signal is As 
and the noise has an amplitude spectral density /-y/rlz, then the signal to noise ratio 
in power, snr is given by 

snr = ^f 1 (6) 
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Figure 4. ASD of noise before line subtraction with the amplitude spectral density of 
the subtracted waveform overlaid - plot (a) . ASD of the noise before (black) and after 
(grey) line subtraction - plot (b). Ratio of the ASD after to before line subtraction 
versus frequency - plot (c). The efficiency of line removal, defined as one minus the 
ratio of the ASD after to before at the line peak, is better than 76% for all the lines. 
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Figure 5. Time series 50 ms of data before (b), and after (d) the subtraction of the 
correction signal (c) generated using the MFT algorithm. So (b) is the sample-by- 
sample sum of (c) and (d). A sine-gaussian having the same sine wave frequency as 
the subtracted line (828 Hz), shown in plot (a), was added to the data in plot (b) to 
show that even bursts at line frequencies are unaffected by line subtraction as long as 
their duration is much less than the response time of the filter. 
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If subtracting the line with a given r leads to addition of noise, the only recourse 
is to increase r, and therefore decrease the bandwidth sampled per Fourier coefficient. 
In the case where the line covers several bins in frequency space, the number of Fourier 
coefficients subtracted will rise as you increase r, resulting in an increase in CPU usage 
proportional to the increase in r. In addition, larger r increases the response time of the 
filter. As discussed in Section 4.6, this will make the subtraction insensitive to a greater 
range of time durations of bursts, which may be desirable from the point of view of signal 
searches. On the other hand, r should not be so large that it exceeds the timescales 
typical of fluctuations in the amplitude or phase of a Fourier component of the line. The 
authors found that r of the order of 8-32 seconds satisfied these requirements for the 
lines detailed in Table 1. Relating signal to noise ratio to efficiency of the line remover, 
lower signal to noise ratio leads to less efficiency in line removal. If the frequency of the 
line is stable and its parameters are slowly varying, one can still increase efficiency by 
increasing r. Otherwise the line subtraction is unavoidably inefficient. However, where 
the main aim of this tool is to suppress lines that make significant contributions to the 
cumulative ASD, and therefore to the time domain amplitude of the data stream, we 
observe that all lines corresponding to discernable steps in the cumulative ASD have 
been efficiently removed after subtraction. 
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Figure 6. Amplitude of sine gaussian signals after they have been injected at random 
positions in a 256 second segment of data, and after the line remover has been run, 
versus the amplitude of the signal before injection. Both amplitudes are expresssed 
in units of the r.m.s. of the data before injections. There were 500 injections with 
amplitudes of between 0.78 and 15.23 times the r.m.s. noise level. As a guide to the 
eye, a straight line of slope one through the origin is also plotted. 
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Figure 7. Response of the line removal algorithm to a large discrete step in amplitude 
of a line during subtraction. The line frequency is 2 Hz. Plot (a) shows the wave fed to 
the line removal algorithm, consisting of the oscillation having an RMS amplitude of 
7mV before 132 s and 14 mV after this time, plus pseudorandom gaussian distributed 
noise with an RMS of 0.2 mV throughout. Plot (b) shows the data after line subtraction 
in the 16 seconds surrounding the discontinuity. The effects of the discontinuity are 
confined to the 8 second interval around the discontinuity. Note that the amplitude 
of the line relative to the noise is far greater than in the engineering data so that the 
effect of the discontinuity is clear above the noise background. 



4-6. Response of the algorithm to injected signals 

Figure 5 shows the response of the subtraction algorithm to a representative sample of 
the test data, but in addition a single short time duration sine Gaussian burst has been 
added to the data. The ratio of the peak signal amplitude to the R.M.S. amplitude of 
the data was 3.125. Plot (a) of this figure shows the added burst. This particular burst 
has been selected as a worst-case scenario, since the frequency of the oscillations inside 
the envelope is 828 Hz, exactly equal to a dominant line frequency in the data. Plots 
(b), (c) and (d) show, respectively, the input data after the test signal has been added, 
the correction signal generated from the data for the same time segment, and the result 
of subtraction of this correction signal. The signal survives in the line subtracted data, 
but the persistent fluctuations from background noise are removed. This separation is 
possible due to the response time of the filter. The Fourier components are measured 
over a time duration of r = 8 seconds, so that any fluctuations in the line parameters on 
a timescale much shorter than r are not significantly corrected by the filter. This is a 
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critical feature of this algorithm for line tracking - the response time of the line tracker 
can be set based on knowledge of the properties of the noise (the stability of the line 
noise background) and the signal, the expected time duration of any transient signals 
that one may wish to avoid suppressing. 

Signal injections were next made at 500 random times in a 256 second data segment, 
where the amplitude of the injection was varied between 0.78 and 15.23 times the r.m.s. 
noise level. The mean separation between injections was approximately 0.5 seconds, an 
average of 16 injections per Fourier transform duration. After signal injection, the line 
subtraction algorithm was applied, then the maximum height of the data in the domain 
of each injected signal was measured. This maximum height was plotted against the 
height of the original injected signal in Figure 6. At small injection amplitudes, the 
maximum height after injection reflects the noise background, since the probability of 
random noise fluctuations under the signal leading to at least one bin that exceeds the 
original injection amplitude increases as the signal to noise ratio decreases. At large 
injection amplitude, the ratio of amplitudes after to before injection fluctuates about 
one. This behaviour is consistent with the maximum height bin in the line subtracted 
data being that having the maximum amplitude in the injected signal, but to this 
maximum height is added a noise term that either increases or decreases the overall 
amplitude by an amount having the probability distribution of the noise background. 

The response of the line removal algorithm to a discontinuous step in the parameters 
of a line was also investigated. Figure 7 shows the results of this investigation. The 
input data, consisting of a 2 Hz sine wave having an RMS amplitude of 7 mV, had its 
rms amplitude stepped discontinously to 14 mV at a time of 132 seconds on this scale. 
Plot (b) shows the output of the filter. Note that the input waveform has been delayed 
by 4 seconds to account for the time delay between the input signal and the correction 
signal discussed in Section 4.2, so that the time axis on both plots represents the time at 
which the signal correction was applied to the data, 4 seconds after the input data that 
generated the correction. The acausality of the filter is clearly visible in plot (b). The 
injection time of the discontinuity was at 128 seconds, corresponding to a correction time 
of 132 seconds, but before this time and after 128 seconds the effects of the discontinuity 
are already apparent in plot (b). Note, however, that the effects of the discontinuity 
are completely confined to the 8 seconds surrounding the correction time, because the 
Fourier coefficients used to evaluate the line parameters are evaluated over r = 8 s. 

5. Summary 

The line subtraction method discussed in this paper successfully removes at least 76% 
of the peak amplitude of all frequency-stable lines to which it has been applied in test 
data. For lines for which a high signal to noise ratio can be achieved, over 99% of the line 
amplitude can be removed. Tests so far have focussed on data acquired at a sampling 
rate of 16,384 Hz, because this is the sampling frequency of fast channels in the GEO 
and LIGO data acquisition systems. At this sampling frequency, the algorithm can 
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subtract 1700 separate frequencies in parallel on a single CPU, and still keep up with 
the incoming data rate. If the rate were 16.384 MHz, at least 1 frequency subtraction 
would still be possible in real time. Because the algorithm operates on each frequency 
component separately, it is amenable to parallelization, in which long as at least 

one frequency can be processed in real time on each CPU, the number of frequencies 
subtracted depends only on the number of CPUs. 

Applications of line-removal in gravitational wave data analysis are wide-ranging. 
This technique could be used as a pre-processing step for existing burst search 
algorithms, and may also be useful in searches for other types of signals. Furthermore, 
it presents the prospect of obtaining a very white spectrum of time domain strain data 
in the gravitational wave channel. If applied successfully, the spectral density would 
be frequency independent. This would allow for more robust time domain searches 
for unmodelled bursts, where the search sensitivity is strongly dependent only on the 
amplitude of a burst signal and only weakly dependent on its shape. The method is also 
applicable to other problems where it is desirable to monitor the amplitude and phase 
of a wave of known frequency, sample by sample, online. Applications in this category 
will be discussed in a future publication. 
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Appendix A. Proof of the EFC evolution formula for Fourier coefficients 

This proof starts with equation 2 for the k th Fourier coefficient, where we separate the 
j = term from the sum and write it at the end. 

F k {xo, xi, ■ ■ ■ ,xn-i) = 2^ Xje N + x . (A.l) 

3=1 

Next write down an expression below for the k th Fourier coefficient of the data forward 
timeshifted by one sample, where the new sample is xn- Rearrange this expression by 
separating the term in xjy and factor out a common phase. 



N 

^ -2,ri(j-l)fc 

Fk{xi, x 2 , • ■ • ,x N ) = l^Xje n 

3=1 

N-1 

■ ^ -2iri(j-l)k -27ri(JV-l)fc 

= x i e N +xn& n 

3=1 



(A.2) 
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The expression in parentheses can be related to ^(xo, xi, ■ • ■ , xjv-i) by using equation 
A.l. 

Fk(xi, xn) = (Fkixo, xat_i) + (x N - x )) . (A. 3) 

This is the desired result. 
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